A light-cone wavefunction approach to open heavy flavor dynamics in QCD matter 
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We calculate the lowest order charm and beauty parton distribution functions in and fragmenta- 
tion functions into D and B mesons using the operator definitions of factorized perturbative QCD. In 
the vacuum, we find the leading corrections that arise from the structure of the final-state hadrons. 
Quark-antiquark potentials extracted from the lattice are employed to demonstrate the existence 
of open heavy flavor bound-state solutions in the QGP in the vicinity of the critical temperature. 
We provide flrst results for the in-medium modification of the heavy quark distribution and de- 
cay probabilities in a co-moving plasma. In an improved perturbative QCD description of heavy 
flavor dynamics in the thermal medium, we combine D and B meson formation and dissociation 
with parton-level charm and beauty quark quenching to obtain predictions for the heavy meson 
and non-photonic electron suppression in Cu-|-Cu and Pb-|-Pb collisions at RHIC and the LHC, 
respectively. 
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I. INTRODUCTION 

The early production of heavy quarks makes them 
some of the most important probes of the quark-gluon 
plasma (QGP) created in ultra-relativistic collisions of 
heavy nuclei [l| . Precise and direct measurements of the 
multiplicities and differential distributions of D and B 
hadrons will soon become available with the vertex de- 
tector upgrades at the Relativistic Heavy Ion Collider 
(RHIC) and at the Large Hadron CoUidcr (LHC). Such 
experimental advances will allow to quantitatively ad- 
dress the key observable in high temperature Quantum 
Chromo-Dynamics (QCD) - the apparent modification of 
particle production by energetic partons traversing a re- 
gion of hot nuclear matter - for heavy quark jets. 
In the framework of perturbative QCD (pQCD), stud- 
ies in this direction have focused so far exclusively, with 
varying degree of sophistication, on the loss of the lead- 
ing parton orparticle energy via radiative and coUisional 
processes 0,0 (commonly known as jet quenching). This 
energy loss leads to a suppression of the observed number 
of energetic particles in heavy ion (A+B) collisions rel- 
ative to their number in proton-proton (p+p, or N-l-N) 
reactions, scaled by the number of binary collisions. 

It is well known that c and h quark partonic level 
energy loss in the QGP is not sufficient to explain the 
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large suppression of non-photonic + e~ measured at 
RHIC Q- These non-photonic electrons predominantly 
arise from the semi-leptonic decays of D and B mesons 
that are usually considered to be formed by heavy quark 
fragmentation outside of the medium. However, due to 
their small formation times, it is possible that D and B 
mesons could form inside the medium and subsequently 
dissociate 0. For this reason, effective energy loss via 
meson dissociation in the QGP 0, H, [§] has also been 
studied in recent publications. Heavy meson dissociation 
approaches have been quite successful phcnomcnologi- 
cally [3], but the question of the possible thermal modi- 
fications to D and B hadrons has remained open. Even 
though such modifications will affect fragmentation only 
if the time scale for their onset is smaller than the in- 
verse virtuality in the parton decay, it is surprising that 
to date there has been no theoretical evaluation of the 
quark distribution (PDF) and fragmentation (FF) func- 
tions at T 7^ 0. 

Heavy flavor provides the ideal testing ground for our 
first quantitative study of such many-body QCD effects 
because the lowest-lying Fock state in the lig ht-cone 
wavefunction expansion for D and B mesons |lOl.[lll| can 
be related to the hadron rest-frame solution of the corre- 
sponding Dirac equation p^.[l3j. Furthermore, potential 
model calculations can be generalized up to temperatures 
of order few Tc [3l using lattice QCD input . In this 
paper we provide the first calculation of the thermal mod- 
ification of PDFs and FFs for open heavy flavor. As we 
will demonstrate below, even in equilibrium with a co- 
moving plasma heavy-light quark bound-state solutions 
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of the Dirac equation persist to temperatures ~ 1.6Tc. 
This provides strong motivation to take the existing stud- 
ies of open heavy flavor dynamics in dense QCD matter a 
step further by combining the partonic and hadronic as- 
pects of the observed high-p^ cross section attenuation. 

For the realistic case of out-of-cquilibrium fast jet prop- 
agation in the thermal medium we combine the partonic 
level quenching with D and B meson dissociation in the 
medium to compute the spectrum of D and B mesons 
prior to their weak decay. We compare the relative im- 
portance of charm and beauty quark energy loss versus 
meson dissociation function of transverse momen- 
tum. We also show results for the non-photonic elec- 
tron suppression in central Au-|-Au and Cu-I-Cu reac- 
tions at y^saT/v = 200 GeV and central Pb-I-Pb collisions 
at y/SNN = 5.5 TeV. We discuss the prospects of testing 
the pQCD theory by upcoming open heavy flavor mea- 
surements at RHIC and the LHC. 

Our paper is organized as follows. In Section |TT] we 
present a detailed baseline calculation of charm and 
beauty quark production, inclusive of the known cold 
nuclear matter effects, and expose the limitations of the 
purely partonic approach to their final-state dynamics in 
the QGP. In Section [1111 we express the parton distribu- 
tion and decay probabilities of heavy quarks in hadrons 
in terms of the light-cone wavefunctions and calculate the 
c and b quark FFs. In Section [TVl we evaluate the heavy 
meson wavefunctions in the vacuum and in the vicinity 
of the QCD phase transition and present numerical re- 
sults for the corresponding distribution and decay prob- 
abilities. A complete calculation that includes D and B 
meson dissociation and c and b quark quenching, relevant 
to current and future heavy flavor/non-photonic electron 
measurements, is given in Section [V] A summary and 
conclusions are presented in Section [VII Background in- 
formation and some of the technical aspects of this work 
are given in the Appendix. Appendix [A] describes nota- 
tion and quantization of parton fields. In Appendix [B] we 
elucidate the relation between the instant and light-cone 
forms of a heavy meson wavefunction. Demonstration of 
the smallncss of the thermal light quark pick-up correc- 
tion to the heavy quark FFs is given in Appendix [Cl Ap- 
pendix |D] presents an analytic solution for the evolution 
of the heavy quark and meson system for time-averaged 
fragmentation and dissociation rates. 



II. A TIMELESS LIMIT 

To correctly assess the discrepancy between the ex- 
pected heavy quark quenching and the measured non- 
photonic electron suppression [l^, improved treat- 
ment of cold nuclear matter effects is necessary. A pre- 
vious study of D and B meson production in p-l-A re- 
actions [l^ indicated that these will affect the nuclear 
modification ratio {Rab) for energetic particles [J. This 



ratio is defined as follows: 
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and is used as one of the main observables to charac- 
terize the properties of the system created in heavy ion 
reactions. In Eq. (|T]), ctab refers to the production cross- 
section in collisions between two heavy ions A and B, 
aN N refers to the cross-section in collisions between two 
nucleons N, and N^^^ refers to the number of binary 
N-t-N collisions in an A-l-B reaction. iVj^^ is determined 
by using a Glauber model [l^. It is, therefore, important 
to establish the role of shadowing. Cronin effect and cold 
nuclear matter energy loss in the analysis of open heavy 
flavor suppression. Note that in our notation p , k • • • 
arc 2D transverse vectors and, for example, px = |p|- 

At the partonic level, the cross sections per elementary 
nucleon-nucleon collision, including the ones for charm 
and beauty quarks, can be calculated as follows: 
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In Eq. ((21) ka , refer to the transverse momenta of the 
two partons a, 6, respectively, that participate in the 
hard collision and Xa , Xh is the fraction of the light cone 
momentum of the nucleon N carried by these. Matted 
denotes the matrix element of the hard collision of a and 
b to give rise to partons c and d. 4>{x, ^) refers to the 
parton distribution functions that give the probability to 
find a parton with light cone momentum fraction x of the 
nucleon light cone momentum. We set the factorization 
and rcnormalization scales to /i/ = /i^ = rriT = {iri^ g + 
py2-ji/2 g^j^j ^j^g Jacobian J reads: 
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subject to the hard scattering constraint kTa,b < Xa,b^fS- 
We have generalized the pQCD coUinear factorization ap- 
proach to heavy flavor production [l^, [2^ to account for 
the non-zero parton transverse momentum distribution 
/(kxa.b) in the hadron wavefunction with {k^^jj)]s[N = 
0.9 GeV^ [2lj. The phenomenological next-to-leading 
order K factor in Eq. cancels in the observable 
Rab{pt), Eq. ([T]). Finally, we note that the momen- 
tum fraction x that appears in Eq. ^ is modified due to 
many-body QCD scattering effects that we discuss below. 

Cronin effect. In p-f-A and A-l-A reactions the Cronin 
effect can be taken into account in the p-differential cross 
sections by including the k broadening of incoming par- 
tons that arises from initial-state scattering [l^iH^. This 
transverse momentum broadening in heavy ion collisions 
is simplest to evaluate if f{kTa.b) is of a normalized Gaus- 
sian form due to the additive dispersion property: The 
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width of the transverse distribution function / grows rel- 
ative to its vahie in N+N coUisions as follows: 
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Specifically, in Eq. © = 0.12 GeV^, Ag = 

{Cf/Ca)\ = 1 fm and ^ ~few is a numerical factor that 
accounts for the enhancement of the broadening coming 
from the power-law tails of the Moliere multiple scattcr- 

Cold nuclear matter- energy loss. Initial-state multiple 
collisions also give rise to radiative energy loss that is 
sensitive to the quark or gluon mass. Theoretical deriva- 
tion of the stopping power of large nuclei for partons 
was given in Ref. [2J| and numerical calculations were 
carried out for the same squared transverse momentum 
transfers and parton mean free paths given above. If fast 

quarks or gluons lose a fraction ea,b = ^e"^'^ '^^ their 
energy prior to the hard scattering, to satisfy the same 
final-state kinematics they must initially carry a larger 
fraction of the colliding hadron momentum and, corre- 
spondingly, a larger value of x. This can be implemented 
in Eq. ^ as follows: 
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in the parton distribution functions 4>a.b/N{xa,bT t^-f)- 
Xa. b will be defined in terms oi Xa ^b in a moment when 
we consider dynamical shadowing. We note that the the- 
oretical calculation of cold nuclear matter energy loss is 
most easily performed in the target nucleus rest frame 
where the incident partons are typically of very high en- 
ergy [13. 

Dynamical shadowing. Last but not least, power- 
suppressed [1^ resummed [1^ [13] coherent final-state 
scattering of the struck partons that leads to shadowing 
in the observed cross sections is included via llSl: 
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thickness of the interaction region is implicit in our ap- 
proach and reflects the actual theoretical findings. 

We now turn to the final-state quark and gluon dy- 
namics that is governed by the formation times of the 
plasma [soj . tq, and of the hadrons, rform 01 ■ Using the 
uncertainty principle, one can evaluate: 
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which grows with underlying jet energy (Et), depends 
on the fractional momentum z = Eh/Eq^g carried by the 
leading hadron and is inversely proportional to the square 
of its mass {nih)- We take the formation time of the QGP 
to be To = 0.6 fm and tq = 0.3 fm at RHIC and LHC, 
respectively, compatible with hydrodynamic studies [2^ . 
Simulations of the early stages of the reaction [l^l suggest 
that the transverse energy density of the partonic system 
builds through a short exponential growth. Therefore, for 
T < To particles will be largely unaffected by the medium. 

For T > To the properties of the soft gluon-dominated 
matter, created in heavy ion reactions with multiplicity 
proportional to the participant density of the colliding 
nuclei: 
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can be constrained by the experimentally measured 
charged particle pseudorapidity density 
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Integrat- 
ing over the position in the transverse plane x in Eq. ^ 
we can determine k from simple parton-hadron duality: 
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Here, dN^ /dy refers to the gluon rapidity density, and 
7/ = — log(tan(6'/2)) is the pseudorapidity, where 9 is the 
angle of the observed particle with respect to the beam 
axis. For example, in central Cu-|-Cu and Au-|-Au colli- 
sions at RHIC, at the average nucleon-nucleon collision 
energy, = 200 GeV, ^ = 250 and 1050, respec- 

tively, and in central Pb-|-Pb collisions at the LHC at 
= 5.5 TeV one expects ^ ~ 2800. The local 
time- and position-dependent temperature for an isen- 
tropic Bjorken expansion can be readily obtained: 
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In Eqs. ([6]) and ^t,u are the relevant Mandelstam vari- 
ables at the partonic level and mc, are the masses of 
the struck target partons. A is the atomic number of 
the nucleus. A consistent picture of cold nuclear mat- 
ter effects requires {i'^)q.,gA^/'^ w {2^i^L/X)q,g in mini- In Eq. glD w e have replaced the proper time T{t,x^) = 
mum bias reactions, which yields ~ 0.12 GcV^ and \/t?' ~ ix^Y ., where x^ is the position coordinate along 



{i^)g « 0.27 GeV2. 

Our results, reported here, bring the treatment of cold 
nuclear matter effects for open heavy flavor production 
on par with their implementation in the study of light 
hadron and direct photon final states psj . Linear depen- 
dence of cold nuclear matter effects on the mean nuclear 



the axis of the collision, by the time t. This is a good 
approximation for x^ small compared to i, meaning ra- 
pidities y small compared to 1. The Debye screening 
scale, ?Ti£) = gT , and the relevant gluon mean free path 
\g = l/crffp with crffs = (9/2)7ra2/TO|,, that enter jet 
quenching calculations can then be evaluated. 



4 



1.2 
1 

^0.8 
< 

cr*0.6 

0.4 
0.2 
0, 



= 200 GeV 



PHENIX R„ 0-10% n , Au+Au 

STAR Rj,p 0-10%/60-80% n+!z'. Au+Au 

Theory, dN^/dy=1050, W-S participant density 

PHENIX R^„ 0-10% Cu+Cu 

Thieory, dN^/dy=250, W-S participant density 







12 



16 



20 



Pt [GeV] 



FIG. 1: Perturbative QCD calculation of the tt" and tt^ + 
7r~ quenching in central Au+Au and Cu+Cu collisions at 
RHIC at -^/snn ~ 200 GeV. Pion suppression data from 
PHENIX m, H and STAR [H is also shown. 



In contrast to the soft gluons discussed above, large 
parton processes follow a binary collision density dis- 
tribution, cx "^-^^^ , that is peaked near the center of the 
collision geometry. In Ref. [3^ a complete and numer- 
ically intensive calculation of jet production and prop- 
agation through the plasma was carried out to confirm 
that the energy loss pattern is determined by the jets 
originating near the peak of the binary collision den- 
sity. Therefore, we here focus only on these hard partons, 
use a running coupling constant at the radiation vertex 
and describe the strength of the interactions of the jet 
with the medium by an effective value ag = 0.28 — 0.32. 
Any uncertainty associated with higher orders in opacity 
and peripheral jets will lead to somewhat larger values 
of as- The energy loss calculation itself is performed 
in the framework of the GLV approach 0, [2J| and the 
Landau-Pomeranchuk-Migdal destructive interference ef- 
fect in QCD is fully accounted for. We evaluate the prob- 
ability distribution P{e) for the parton to lose a fraction 
of its energy loss e = due to multiple gluon emis- 

sion. The quenched quark or gluon spectrum is then 
readily obtained: 
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In Eq. the factor (1 — e)^^ arises from the transverse 
phase space Jacobian Id^p/d^pQuenc/ij^ Hadronization is 
performed using the appropriate vacuum fragmentation 
bmction Dii/q g(z, fifr) as follows: 
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FIG. 2: Perturbative QCD calculation of the suppressed D 
and B meson production cross sections in central Au+Au (top 
panel) and Cu+Cu (bottom panel) collisions is compared to 
the magnitude of light pion quenching at y^sjvjv = 200 GeV. 
The effect of parton mass on the cold nuclear matter effects 
is illustrated via uncertainty bands. Non-photonic electron 
data from PHENIX ^ and STAR ^ is also shown. 



where z refers to the fraction of the light cone energy of 
the initial parton q, g, carried by the hadron h. 

Eqs. ([2]), (|12p and (jl3p summarize the traditional ap- 
proach to jet quenching that has been very successful 
in the description of the observed attenuation of light 
hadron production in A+A reactions. Figure [1] illustrates 
the level of agreement that can be achieved between 
the perturbative QCD theory, described in this section, 
and the experimental measurements of [ssl . [3^ and 
i(7r+ + 7r-) [Hi at RHIC. In central Au+Au and Cu+Cu 
reactions we used gluon rapidity densities dN^/dy = 
1050 and dN^/dy = 250. Identical calculation for 
open heavy flavor production in heavy ion collisions at 
y/SNN = 200 GeV is presented in Fig. O It is easy to 
see that Cronin enhancement plays an important role at 
intermediate px ^ 4: GeV. The bands show the sensitiv- 
ity of cold nuclear matter effects when we vary all parton 
masses in the range from to mc,mb, respectively. At 
mid-rapidity, y = 0, the effect of high-twist shadowing is 
small and the uncertainty band is dominated by initial 
state energy loss. Cold nuclear matter effects amplify the 
disparity between the heavy and light parton quenching. 
We have schematically illustrated the latter by showing 
the magnitude of pion attenuation from Fig. [1] 
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In contrast, the experimental results on non-photonic 

electrons [IB, [13 R'aa ~ ^^A^ Pt > 5 GeV are in clear 
contradiction with the small quenching of B mesons that 
give an increasingly important > 50% contribution to 
non-photonic e+ -|- in this region 0, H^. An insight 
into the possible cause for this well-known by now dis- 
crepancy can be gained by examining Eq. Q . The small 
pion mass ensures that the parent light quarks and glu- 
ons fragment outside of the QGP in accord with the tra- 
ditional picture of jet quenching. The large D and B 
meson mass, however, implies that charm and beauty 
quarks will fragment inside the hot and dense medium: 
Tform oc Consequently, the competition between 

heavy meson dissociation and the c and b quark decay in 
the QGP is a likely physics mechanism that may natu- 
rally lead to attenuation of the beauty cross section as 
large as that for charm 



III. DISTRIBUTION AND FRAGMENTATION 
FUNCTIONS OF HEAVY QUARKS 

Simulations of heavy quark fragmentation and dissoci- 
ation in the QGP require knowledge of the correspond- 
ing parton decay and distribution probabilities. These 
are related to the wavefunctions of their parent or decay 
hadrons that can be represented as follows: 

oo n 

|P+;J)=at(P+;j)|0)= ^ In 

n=2(3) 
/ n \ 

X ilj{xi, ki; a,) X! 

Vj=i / Vi=i 

X I • • • i^q, + kg, ) "9. ) • • • (Xqj + kg, , Qfg J ' ' ' 

(14) 




Here, P+ = (P+,P) are the large light-cone mo- 
mentum and transverse momentum components of the 
hadron, and the momenta of the partons arc given by 
{xiP^ ^XiP 4- ki). Our conventions for the quantization 
of the fermion and vector fields, including notation and 
commutation relations for the creation and annihilation 
operators a^, a, • • • , are given in Appendix El In Eq. (jl4[) 
ttg; is a set of additional relevant quantum numbers, such 
as helicity and color. From the normalization of the me- 
son state of fixed projection A of the total angular mo- 
mentum J: 

(P+; J\P+'- J) = 2P+{2'Kf5{P+ - P+')5'(P - PO-^AA' , 

(15) 

an integral constraint on the norm of the light-cone wave- 
functions ip{xi^ ki; ai) can be derived. 

In Eq. (|14[) proper construction of the non-perturbative 
lowest order Fock component in the expansion of a 
hadronic state in a quark and gluon basis, qq for mesons 
and qqq for baryons, appears to be the most impor- 



tant. When properly constructed, higher order compo- 
nents that are generated via parton splitting, will have 
the correct quantum numbers through preservation of all 
relevant global symmetries in QCD. These components 
become dominant in the small momentum fraction x re- 
gion of sea quarks and gluons. On the other hand, for 
X ~ C(l) the momentum distribution of valence quarks 
can be determined from the rclativistic quark model. For 
heavy quarks (Q), pcrturbative splitting is suppressed 
by factors ©(Ag^^/Afg) relative to hght quarks [q) and 
gluons {g) and if we focus on the charm and beauty dis- 
tribution and fragmentation functions the lowest order 
Qq (or Qq) Fock component for heavy mesons will be a 
good starting point for our calculation: 



|P+;J) = 
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{xP+ + k)6t _ ^)p+ _ k) |0) . (16) 



In Eq. si, S2 refer to the spins of Q and q and we 
have used the notation: 



1 

-1 



(17) 
(18) 



to ensure the correct spin structure of pseudoscalar ^5*0 
and vector ^Si states and proper normalization [37| . A 
set of polarization vectors for J = 1 in Eq. (fT5)) in the 
meson rest frame is given by: 




,£3 



(19) 



The color singlet structure is ensured by 5ciC2 where 
ci, C2 are the color indices of Q and g, respectively. Fol- 
lowing Eq. (jlSp we can derive the normalization for each 
helicity A: 



2(27r) 



dxd^k |V^(x,k)|2 



1 . 



(20) 



Note that the meson wavcfunction can also be written in 
an instant form. For example, in the rest frame of the 
meson: 

where q represents the relative momentum between 
quarks and f{q) is the Fourier transform of the meson 



wavcfunction in position space. E}^ = y P^ + nif^ is the 

meson energy and Ei = \/ qf + m? are the quark en- 
ergies. Conversion between the instant and light-cone 
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forms is facilitated if harmonic oscillator, or Gaussian, 
approximation is made for the shape of the wavefunction 
f{q) [sH], which allows to easily ensure identical trans- 
verse width (k^) 01 for both forms. 

In the light-cone gauge = 0, the heavy quark dis- 
tribution function is given by [s^ : 



dy 
27r 



-ixP y 



x(P+|^Q(y-,O)^^Q(0,O)|P+) , (22) 

where ij^q is the heavy quark field operator, IP"*") are ap- 
propriately normalized hadronic states with large light- 
cone momentum P"*", and the average over the hadron 
spin/polarization states is implicit. For both the pseu- 
doscalar and the vector cases, substituting Eq. ([TO]) in the 
perturbative QCD definition, we proceed to a straight- 
forward evaluation at tree level. We first carry out the 
contractions of the heavy and light quark creation and 
annihilation operators, perform the color sums and the 
integration over the internal kinematic variables to arrive 
at an intermediate result: 



d^k dxQ f dy 



-i{x-XQ)P+y 



(27r)3 2xQ J 2n 
xK'(xQ,k)p^^^(A/(A)A/t(A))^,^^ 



where (p) are Dirac spinors defined in Appendix [X] 
With the polarization average yielding (5siS2 ^^id the spin 
sum resulting in xqP^ we readily obtain the final result: 



(I^Q/hix) 



2(27r) 



dxQd^k\^P{xQ,k)\^dixQ-x) . (23) 



In contrast, the tree- level calculation of fra gme ntation 
of the heavy quark into hadron h, defined as (39[: 



„ , ^ f dy ,p+ - 1 1 7+ 

Dh/Q{z) = z -^e ' " -Trco/or^TrDirac-^ 

x{0\^{y-,O)aliP+)ah{P+)mom , (24) 



can be shown to be proportional to (0|5t (• • • )feg( - • • )|0) = 
and, hence, vanishes. This result is also reinforced 
by the mismatch at tree-level between the large light- 
cone momentum of the parent heavy quark and the cor- 
responding momentum of the heavy quark inside the me- 
son. Since both z = /pg < 1 and < 1 , in the region 
of interest D/i/q(0 < z < 1) = 0. 

The first non-trivial contribution to Di^^q{z) comes 
from the following matrix element that replaces the tree- 
level expression in Eq. 



(0|Vq(2/-,O) J d%{-igA'^^^QrT''i'Q) j d'y2{~igAZ%i^T'^i^,) x a^Q(PimP'i)W2)aQ{P2) 
X / dS:i{-igAli;,^^T^%) [ d^yi{-zgAli,Qj^Tf'^Q)^{0,0)\0) , 



(25) 



and corresponds to the process with Feynman diagram 
shown in Fig. [S] T" are the standard GcU-Mann ma- 
trices. Our approach is similar to the one outlined 
in [33, [13, El but in the calculation we retain q, the Qq 
relative momentum distribution, in the matrix clement. 
This will allow us to calculate the lowest order correction 
in q to Dh/Q[z). 




After converting the integral over the transverse mo- 
mentum of the outgoing free light parton pq into an inte- 
gral over the virtuality of the heavy quark we obtain the 



FIG. 3: Perturbative QCD calculation of heavy quark frag- 
mentation matched to the lowest lying Fock state in the heavy 
meson wavefunction. 
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heavy quark fragmentation functions: 



Dh/Q{z) 



(27r)32^xi(l-xi) (2^)32^X2(1 -X2) 



M(j)s,.; M(j).s. 



Tr 



V2 



7 



V2 



1 - z 



1 -Pa-mq 

Vs'Sp'2)Us2{P2)1^ 

1 



1-Pb-m, 



TV[7+(7.p)] 



(26) 



Here, Ilf^^{pa - Pq) = i(-.9pi' + {nfiK + n^k^)/{n ■ k)) 
is the gluon propagator, with = [0, 1, 0] in Ught cone 
coordinates, and s — ~ p^, with Pa — Pb = Pq- If we 
define: 



T(z) = M(j),,,,M(j)s 



ds9 s-^ 



1 - z 



Tr 



7 



7 • Pa - "IQ 
Vs!,{p'2)Us2 {P2h^ 



7^'"si (pl)«6-; (^'1)7'' (7 • + "^9)7° 



1 -Pb-mq 



-n„.n, 



Tr[7+(7-p)] ' 



(27) 



the heavy quark fragmentation functions into heavy 
mesons are then given by: 

D.M^)=l'\¥^'^^i^nz). (28) 



(27r)3 



(27r)3 



Note that the connection between the instant-form wave- 
function f{qi) and the hght-cone wavefunction '4'{xi, ki) 
is discussed in detail in Appendix IB] Thus, Eq. ([^5]) is 
equivalent to Eq. ((26|) . 

The spin sums for the Qq (P) and 3^1 (V) states 
can be derived, yielding: 

UsApi)^s'Sp'iW{P)s,s', 

^ (:^i+toq)(1+7°)(/i +™g)7^ 
v/2(ri' + '"Q)V2(pf + m,) ' 

^ ".i(Pi)«4(p'i)M(^).i4 

^ (j^i+mQ)(l+7°)(/i+m,)^p(A) 
v/2(p'j' + mQ)v/2(pf + m,) 



(29) 



(30) 



where ep(A) is the 4-component polarization vector for a 
specific helicity projection A of the ^Si state. In deducing 
the above expressions, we can first compute the spin sums 
in the rest frame of the heavy meson, then generalize 
them to an arbitrary frame. 

Let us first neglect the effects of quark motion on the 
meson mass and assume in the rest frame that M = 



mq + rUq, r = niq/M, 1 — r = ttiq/M. In a Lorentz 
boosted frame the light and heavy quark momenta can 
be written as: 



Pi = (1 - r)p + qip , p\ ^rp- qip . 



(31) 



In Eq. ()3ip qip is the boosted relative momentum be- 
tween the quarks in the lowest Fock state {q at rest). A 
classic approximation [s^, SOj EH is to assume that the 
wavefunction f{q) of the heavy meson as a bound state 
of two quarks goes rapidly to zero when q ^ Q. The dom- 
inant contribution is, hence, given by the region g ~ 
and in this limit the spin sums in Eq. (|29p and Eq. (^0]) 
can be reduced to: 



UsAPl)^s',{p'l)M{P)s,s', 



^«.,(pi)«,.(p'i)A/(F),,,. =- 



M 
M 



V(^-M), 

(^ + A/)/p(A), 
(32) 



which are the same as the ones used in Refs. [stI I40ll4l|. 

In this paper we employ the complete expressions for 
the spin sums, Eqs. ((29l) and ((30l) . and calculate T{z) de- 
fined in Eq. (|27p . This evaluation is technically very com- 
plex, with C(10*) terms, and requires the use of packages 
for Dirac algebra manipulation 44] . Having obtained the 
full result, we can take the g — > limit to obtain the low- 
est 0{q^) contribution to heavy quark fragmentation into 
the ^Sq state: 



(z) 



(Z- l)2z 



3(r-l)2(2(r-l)r + l)2^ 



3M3r2((r- l)z-M)6 
+2(r - l)(r(18r - 19) + 6)z^ + (GSr^ - 74r + 21) z^ 



+ {36r-18)z + 6 

and the one into the 35*1 state: 
{z^lfz 



(33) 



Tl"\z) 



M3r2((r- l)z+ 1)6 
z ((r - l)^(2(r - l)r + 3)z^ + 2{r ~ l)(r(2r - 1) -f 4)^^ 

(34) 



+6r(2r - l)z + 9z + 4r - 6) + 2 ^ , 



respectively, which recover the results given by [33, |41| . 

Next, we proceed to evaluate the 0{q^) correction to 
the fragmentation function. To do so, we not only select 
the terms with 1-st power of qip and q2p, but also take 
the large pg limit (keeping the leading power of Pg)- We 
are able for the first time to obtain a meson structure- 
dependent correction to the heavy quark fragmentation 
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functions. For the ^Sq state or result reads 



q+P q+p AP 



+32]z2 + ((25 - 8r)r - 15)z + 3 [> , 
while for the ^Si state we get: 



(Z- 1)2^2 



{( 



6M3(l-r)r3((r-l)z + l)8 

^ (122^ + 1922 + 62 -16)' 



82'' (22 + 2 + 2)^ 



+22 (182^ + 502^ - 6822 + 5I2 - 36) 



-1762^ 



20322 - I2I2 + 16) r2 + (2 - 1)2 (32^ + 492^ 



^592- 



+ 5O2 - 6) r - (2 - if (II22 - 132 + 9) ) 



M2 



3rS2S + rM-622 + 132 



9iP 

14) z^-r^ (122^ + 272^ 



+292 - 32) 2^ +r^ (482" - 162^ + 2O22 - 92 - 28) 22 

-r2 (572^ - 8O2* + 262^ + 5722 _ ^ -^g-j ^ 

-(2 - i f (62^ - 22 - 42 + 6) + r{z - 1)2 (3O2'' - 92^ 



-IO22 + 322 




(36) 



In Eq. ([M)) we observe that for the triplet state the co- 
efficients of qip and q2p are different, which results from 
the different combinations of 7 matrices. For example, a 
typical term in the intermediate results for the ^Si state 
may look like 

Tr[A^ip(l+7")(^^ + il/)B(:^ + M)] , 

where A and B are polynomials of Dirac matrices. This 
is not always equal to: 

Tr[A(:^ + A/)B(|^ + M)(l + 7")^2p] , 

even after exchanging the momenta (1 <-> 2). On the 
other hand, qix, q^p (or qui the rest frame) are them- 
selves internal integration variables with equal weights 
'0(a;i,ki) (or f{qi) in instant form), see Eqs. (|26p and 
(PS)) . Hence, Eq. ([55)1 can be written in a symmetric 
form if desired. 

Finally, we note that if the relative momentum of the 
light and heavy quarks is not neglected then in a boosted 
frame the large momentum of the meson in Eq. (j3ip 
will be partitioned proportional to the parton transverse 



3(r - l)^r2^ + (r - lf{2r + 5)2^ + 3(r - lf{Ar - 7)2^ 
^(r - l)(r(r(8r - 63) + 88) - 36)2^ + [r((81 - 28r)r - 90) 



(35) 




FIG. 4: Fragmentation functions of charm and beauty quarks 
into and Si heavy meson states. Solid lines are for negli- 
gible momentum spread between the parton constituents of D 
and B. Dashed lines represent the first calculation for if 7^ 0. 
Insert shows the fractional correction that arises from the 
structure of the final-state hadrons for both the pseudoscalar 
and vector states. 



mass, rriTQ = \jf^Q + (1^2), mT q ~ y^'^q + ^) ■ In this 
more general case M = niTQ + rtiTq and r = rriTq/M, 
1 — r = ttitq/M. In this work we calculate (k2) = (q2) 
from the light-cone wavefunctions of heavy mesons and 
use the corresponding values of r both in the vacuum 
and at T 7^ 0. This allows for a more consistent descrip- 
tion of heavy quark fragmentation without treating r as 
a phenomenological parameter. 

With expressions for T(z) at hand, after integrating 
over the relative momenta of Qq system in the rest 
frame of the heavy meson, we obtain the 2-dependence 
of heavy quark fragmentation functions for the ^Sq and 
^Si states. Keeping in mind that the relative momen- 
tum is rather small and in the calculation |g| <C rM 
and \q\ <C (1 — r)M was implicit, we choose the upper 
bound of this integration to be rM/2 due to the fact that 
r < (1 — r). The absolute normalization of Di^^q(z) has 
been fixed by the known branching ratios of the charm 
and beauty quarks into heavy mesons and baryons, such 
that J2iIo Dh,/Qiz)dz = 1. 

Numerical results for c and b quark fragmentation into 
the ^So and the '^Si D and B mesons states, respec- 
tively, are shown in Fig. U) The solid lines in the figure 
denote the standard 0{q^) calculation of Dii/q{z), in ac- 
cord with [13, Ell . The dashed lines signify the addition 
of the 0{q^) correction, the first correction that arises 
from the internal structure of the final-state meson. The 
magnitude of the 0{q^) term is best seen in the insert 
of Fig. H in the ratio D^°+^^z) / D° (z) . While on an ab- 
solute scale the correction to the shape of heavy quark 
fragmentation can be ±50%, it doesn't change apprecia- 
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bly the position of the large-z peak and, hence, will not 
alter the cross section for the observed D and B meson 
production. To facilitate the comparison between the 
FFs in the vacuum and at T 7^ we will only use the 
leading result for structureless hadrons, Eqs. p3)) and 
(131. 



THE FATE OF HEAVY MESONS BELOW 
AND ABOVE Tc 



IV. 



In this section we move beyond the vacuum T = 
case. Potential models have been quite useful in study- 
ing the non-perturbative aspects of QCD. In particular, 
for quarkonia one can put the potential model on firm 
grounds by using an effective field theory (NRQCD) [43| . 
More generally, heavy-light mesons like the D and the B 
and their excited states at zero temperature have been 
studied in detail [l^, by treating the light g as a Dirac 
particle moving in a confining potential set up by the 
heavy quark Q. In this section we investigate the bound- 
state wavefunctions for D- and B-mesons at rest and in 
equilibrium with a thermal medium or in a co-moving 
plasma. These can be employed to calculate the PDFs 
and FFs for quarks/mesons, as discussed in Appendix IB] 

Lattice QCD results for the free energy of a system 
at temperature T containing two infinitely heavy quarks 
separated by a distance r [1^ have been used to extract 
the static potential between these quarks: 

r -^ + <Jr, r<r,ned{T) 

(T) , r > 

' mea 

(T) 



V{r) 



ai{T) cxp(-M(T)r) 



(37) 

Here, rmed{T) signifies the temperature-dependent dis- 
tance scale at which the medium breaks off the linear con- 
fining potential, ai{T) represents the effective coupling 
at a given temperature, and /^(T) is the Debyc screen- 
ing mass. Interpolation is typically used to smooth out 
the sharpness in the potential at r = ^med and we de- 
scribe the interpolation below. This potential has been 
employed previously to study the fate of heavy-heavy 
bound states in the QGP formed at temperatures above 
Te [illil. 

For future convenience, it is useful to separate the "vec- 
tor" {Vv{r)) and "scalar" {Vs{r)) parts of the potential in 
Eq. (I37D. This is done as described for T = in [H, [11] . 
In Eq. p7[) the term proportional to cr, which in the ab- 
sence of a medium would be a confining linear potential, 
is taken to be scalar. The Coulomb part is naturally in- 
terpreted as the zeroth component of a vector potential. 
Thus we consider: 



^vij") { Qi cxp( — /jr) 



r < r. 



Vs{r) = 1 



r > r, 
d 



med 
7ned 



(38) 



a 



The typical parameters used in Eq. ([37]) are 
= 0.445, ai(r) = 375'(27rT), and fi{T) = 



1.417y^l + ^g{2TTT)T. For T > l.lTc, we use the per- 
turbative two loop expression for g(27rT). For T be- 
low O.OTc, where the perturbative expression for g{2'nT) 
does not work, we take g — 2.05 which corresponds to 
a = 0.445; and near Tc we interpolate between the two. 
The strength of the confining potential is taken to be 
(7 = 0.224 GeV^ and = 192 MeV. The value of a is 
w 5% larger than in [3] and we take Mc = 1.34 GeV and 
A/fc = 4.79 GcV. r^ned = min(0.4^ 11) fm, where 1.1 fm 
is the typical string breaking scale [4J] . Vs (r) is smoothed 
out on a scale much smaller than rmed{T)- Vy{r), follow- 
ing [il] , is smoothed by interpolating between a/r at r = 
rmed(T) and ai exp(— /zr)/r at r = ti{T), where ri(T) 
is larger than rmed{T). For our calculations, we take 
ri{T) = 2r™ed(r) for T < l.lTc and ri(T) = 1.25 fm for 
T > l.lTc. This gives a slightly stronger potential than 
that used by . 

The presence of bound states of heavy-heavy sys- 
tems, like charmonia and bottomonia, can be explored 
by solving the non-rclativistic Schrodingcr equation for 
the heavy quarks in the presence of the potential Eq. (|37p . 
Here, we only look at the / = wavefunction, for which 
the equation has a form, 



1 1 d^r^ 
2{mQ/2)r dr^ 



V(r)V' = E-iP{r) 



(39) 



where mg is the mass of the heavy quark. We find that 
J/ip bound state solutions exist up to about 2Tc though 
they have a binding energy, quite sensitive to the pre- 
cise value of a, much less than the temperature. The T 
disappears at much higher temperatures but its binding 
energy becomes smaller than T near 1.5Tc. These values 
are larger than the values reported in [4^ because of two 
reasons. First, the value of cr is a little larger. Second, 
we use a slightly different interpolation to smooth the 
potential between r < r,„ed and r > r„ied- Table [T] shows 
the binding energy and sizes for the J/iJj and the T at 
various temperatures above Tc. 

To investigate the existence of heavy-light mesons in 
the QGP we need to solve the Dirac equation for a light 





T (GeV) 


Et (GeV) 


(r) (fm) 


1.2 


0.230 


0.042 


0.722 


1.4 


0.269 


0.030 


0.804 


1.6 


0.307 


0.024 


0.869 


1.8 


0.346 


0.020 


0.923 


2.0 


0.384 


0.017 


0.968 


1.2 


0.230 


0.344 


0.222 


1.4 


0.269 


0.301 


0.232 


1.6 


0.307 


0.273 


0.241 


1.8 


0.346 


0.254 


0.248 


2.0 


0.384 


0.241 


0.254 



TABLE I: The binding energy and mean radius for J/i/; (up- 
per table) and T (lower table) as a function of the tempera- 
ture. 
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particle moving in the potential created by a heavy quark: 

H'il) = E^ , (40) 

where ^ is the four-component Dirac spinor and H is the 
Dirac Hamiltonian. H has the following form: 

H = -if -I'd, + V„{r) + 7"(7n + Vs{r)) . (41) 

The eigenfunctions of the spherically symmetric Dirac 
Hamiltonian can be written as [H, [ij, 0, : 

where J^jf^ are angular functions that are obtained by the 
angular momentum addition of s = i to Z giving j . The 
energy eigenequations for the radial wavefunctions F(r) 
and G{r) that we need to solve are as follows: 

F'{r) - !^F{r) = {-E + V^{r) + m + Vs{r))G{r) 

G'{r) + ^G(r) = {E - V.,{r) +m + Vs{r))F{r) , 

(43) 

where k = ±(j + 1/2). For the lowest energy states that 
we are interested in, j = 1/2 and k = —1. We solve 
Eq. (|43|) numerically using the "shooting method" [13]. 
Since the potential tends to a constant at large distances, 
we begin from an exponentially decaying solution at a 
very large distance. We solve the differential equation 
implied by the Dirac equation for this boundary condi- 
tion and choose the value of the energy E such that the 
value of the wavefunctions F and G at ?■ = is 0. 

Since the heavy-light mesons are larger in size when 
compared to quarkonia, it is natural to expect that they 
will be affected more severely by color screening. Indeed, 
for a light current quark mass of about 0.005 GeV, we 
don't find even very weakly bound states for T > Tc- 
However, this conclusion is quite sensitive to the param- 
eters chosen above. For example, by simply changing 
rmed from 0.4^ fm to 0.45^ fm [Ij] we find bound 
states up to ~ 1.2Tc, albeit with binding energies Et 
about 50 times smaller than T. A quark moving in the 



T (T,) 


T (GeV) 


Et (GeV) 


y(r2) (fm) 








0.730 


0.468 


0.2Tc 


0.038 


0.733 


0.466 


OAT, 


0.077 


0.611 


0.464 


0.6Tc 


0.115 


0.256 


0.501 


0.8Tc 


0.154 


0.098 


0.632 


LOT, 


0.211 


0.043 


0.785 


1.2Tc 


0.230 


0.031 


0.970 


1.4Tc 


0.269 


0.017 


1.263 


1.6Tc 


0.307 


0.009 


1.636 



TABLE II: Properties of the D and B meson bound-state 
solutions taking the effective light quark mass to be m(T) /\/2- 
These persist up to temperatures ~ l.BTc. 



T (Tc) 


1 (OeV j 




Vv / (™) 


1 HT 


911 


u.uou 


U.D 1 u 




0.230 


0.047 


0.797 


l.lTc 


0.269 


0.032 


0.933 




0.307 


0.024 


1.044 




0.346 


0.020 


1.119 




0.384 


0.018 


1.158 



TABLE III: Properties of the bound D and B meson states 
taking the effective light quark mass to be m{T). 



QGP will have an effective thermal mass larger than its 
bare mass niq. For T > Tc perturbative calculations 

give m(T) = y%[5(27rT)r]%m2. For T < T^ this 
form is not valid and we simply interpolate the pertur- 
bative value above T^ to the bare mass at T = 0. One 
would expect the effective mass of a quark moving un- 
der the influence of a potential in the thermal medium 
to lie somewhere between niq and m[T). For example, 
results obtained by solving the Dirac equation for the 
light quark mass taken to be m(T)/V2 are given in Ta- 
ble [TTl We remark that the very small rise in the binding 
energy that can be seen in Table |TT] as we increase the 
temperature from T = to 0.2rc is an artifact of the 
interpolation we have used for the thermal mass below 
Tc and is not numerically significant to affect our results. 
The use of light quark mass equal to m[T) gives bound- 
state solutions with weak binding energies up to even 
above 2Tc, see Table IhH We also establish the existence 
of D and B meson bound-state solutions well above the 
phase transition temperature for constituent light quark 
masses and/or for the "maximum confining potential" 
compatible with lattice data: all the way up to ~ 2Tc for 
ruq = to(T)/-\/2 and ^ 3Tc for ruq = m{T). 

The bound-state wavefunctions, calculated at temper- 
atures T = OTc, 0.4Tc, 0.8Tc, 1.2Tc with the light quark 
mass taken to be m{T)/^/2 are employed to evaluate the 
PDFs and FFs of heavy quarks in a co-moving thermal 
medium from Eqs. (|23p and (|26|) and our results are pre- 
sented in Fig.[5l The distribution functions 0q///(x) be- 
come narrower in x as we increase the temperature, which 
is intuitively expected because the wavefunctions become 
broader in position space. In contrast the peak position 
and width of the decay probabilities are determined by 
the boost parameter r. These remain nearly constant 
because the decrease in the mean transverse momentum 
squared is largely compensated by the growth in the ther- 
mal mass. Hence, we have used an arbitrary scale factor 
s to better separate the curves in the bottom panel of 
Fig. O We remark that the double-peak structure arises 
from the decay of the excited ^Si states into ^5*0 states 
which effectively shifts the momentum fraction z to lower 
values. 

To summarize, as a rule we find that D and B me- 
son bound-state solutions persist above Tc- Their small 
binding energy and large radius, however, will greatly 
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20 



m^' = m'(T)/2, T<1.5T^ 




X, Z 



FIG. 5: The evolution of PDFs (pQ/ni^) (top panel) and the 
FFs Dh/q{z) (bottom panel) with T. Blue/black lines refer 
to charm while red/orange lines refer to beauty. A scale factor 
s has been used to better separate the decay probabilities. 



facilitate their subsequent dissolution in the presence of 
interactions. 



V. APPLICATION TO HEAVY MESON 
PRODUCTION IN HEAVY ION COLLISIONS 

Heavy flavor dynamics in dense QCD matter critically 
depends on the time scales involved in the underlying 
reaction. Two of these timescales, tq and Lqgp, can be 
related to the nuclear geometry, the QGP expansion, and 
the properties of bulk particle production [3^ . They sig- 
nify the onset and disappearance of QGP-induced effects 
and determine, in relation to the formation time Eq. ([5]), 
whether the suppression of the observed cross section 
arises predominantly at the hadronic (rform ^ Lqgp) 
or partonic (rform > Lqgp) level. As discussed in the 
previous section, the dissociation of D and B mesons in 
the vicinity of Tc can be facilitated by their small binding 
energy, or, equivalently, their broad wavefunction in co- 
ordinate space. Whether such thermal effects take place 
in practice, however, depends on the time they need to 
develop. We can roughly estimate thi s tim e by boosting 

from TableHH 



the expanded size of the hadron, 

by the 7 = 1/Vl — v'^ factor. When compared to rform, 
which is determined by the virtuality in the parton decay 
process, this time is large and suggests that the fragmen- 
tation component of the heavy meson dynamics in heavy 
ion collisions may not be affected by the QGP. In what 
follows we will study this "instant wavefunction limit" in 



detail. Before we proceed we remark that in Appendix [Cl 
we have also estimated the correction to the FF D}^/q{z) 
from the pick-up of a light thermal quark strictly in the 
fragmentation region < z < 1 and found this effect 
to be small, < 3% and 15% for D and B mesons at 
Pt > 4 GeV, respectively. 

A meson that is formed and propagates inside the 
medium will undergo coUisional broadening and disso- 
ciation 0, Q . The dissociation rate can be evaluated as 
follows: 



<,rdi; 



s(pT,t)) 



d_ 

dt 



ln(l - P,(pT, mg, t)) , (44) 



where the survival probability is given by the overlap of 
the initial wavefunction t/jo with the wavefunction at time 



P. 



"XT 



S'k.dx V'}(k, a;)V'o(k, x) 



(45) 



Details of the diagrammatic calculation and resummation 
of the interactions that lead to the final-state broaden- 
ing of the D and B meson wavefunctions, ipoiya^x) — > 
ip f (k, x) , can be found in (3] . We here remark that the 
harmonic oscillator wavefunction form, discussed in Ap- 
pendix |b1 allows for an analytic evaluation of this broad- 
ening. As in Eq. ^j-^^ is proportional to the cumu- 
lative transverse momentum transfer from the medium 
to the propagating system and is here calculated for 
Bjorken-expanding QGP. Formation rates are calculated 
by taking into account all final-state heavy meson states 



1 



(rform (PT,0) 



V / dzDH,/Q{z) 
, Jo 



'T-formiz, PT,mQ,t) 



(46) 



Employing the above calculated rates, Eqs. (|44|l and 
(|46p . the concurrent processes of c and b quark fragmen- 
tation and D and B meson dissociation are described by 
the following set of rate equations 0: 



dtf'^iPT^t) = - 



1 



rform (pT,t)) 



f'^iPT^t) 



{rM^,t))l d^^<^QM-)f"iPTMA^7) 



dtf"{pT,t) 



{Tdiss{PT,t)) 



f"{pT,t) 



dz-DH/Q{z)f'^{pT/z,t) .{A8) 



{T{ormipT/z,t)) Jq Z 

Here, f^^ipT^t) {f^{pT,t)) is the differential cross sec- 
tion to find a quark Q (hadron H) with a fixed rapidity y 
and transverse momentum pr, at a time t. We have sup- 
pressed the rapidity dependence for clarity of notation. 
Note that the reason for which the asymptotic t —f 00 
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DC 



< 
< 
DC 



1 - 

0.8 - 
0.6 - 
0.4 - 
0.2 - 

— 

1 

0.8 - 
0.6 - 
0.4 - 
0.2 - 



s -200 GeV 



Central Au+Au 



5 = 3-5 



■ — B mesons 
— D mesons 




Raa light 



Instant wavefunction form 



—I 1 

s"^ = 200 GeV 



Central Cu+Cu 




Instant wavefunction form 

J I I , 



10 15 
Pt [GeV] 



20 



FIG. 6: Suppression of D and B hadron production from 
meson dissociation and heavy quark quenching in central 
Au+Au (top panel) and Cu+Cu collisions (bottom panel) at 
= 200 GeV at RHIC. 



solution will exhibit suppression of the cross sections is 
that both fragmentation and dissociation processes emu- 
late energy loss by shifting the quarks/hadrons to lower 
momenta. Under certain simplifying assumptions the 
system of equations, Eqs. ((47|) and (|48|) . can be solved 
analytically, as shown in Appendix [Dl 

At present, there is no reliable way of incorporating 
the fluctuations in partonic energy loss in rate or trans- 
port equations. Therefore, we include the early-time 
heavy quark inelastic scattering effects approximately as 
a quenched initial condition: 



f"{PT.t) 



d 



/^(PT,0) = 0. 



dyd'-^p 



dyd'^p 



(49) 



Here, the attenuated partonic spectrum "^"^^^ip — is cal- 
culated differentially versus pt using Eq. (fT2)) . The rel- 
evant mean quenching time - the time that the physical 
system of interest spends in a quark state - can be esti- 
mated from the analytic solution for /'^{pr, t), f^{pT, t) 
given in Appendix |D] On the other hand, Tform repre- 
sents the important stage until the first fragmentation 
of the heavy quark into heavy meson. We use the aver- 
age of these two results to calculate the initial condition, 
Eq. (091). 

We integrate numerically the above set of coupled or- 
dinary differential equations and use the same initial soft 
gluon rapidity density dN^ /dy as in the simulations of 
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FIG. 7: Nuclear modification for the single non-photonic elec- 
trons in central Au+Au and Cu+Cu collisions at RHIC. Data 
is from PHENIX [11 and STAR [Tl coUaborations. 



tt" quenching. The corresponding suppression of open 
heavy fiavor in ^'snn = 200 GeV central Au+Au and 
Cu+Cu collisions at RHIC is shown in the top and bot- 
tom panels of Fig. [6j respectively. For D mesons the 
Cronin effect is clearly visible around pT ~ 4 GeV and 
this is a notable difference from our previous study 
where initial-state kx diffusion was not included. For B 
mesons the change in heavy quark velocity /3q, which 

_ 2 

controls the Qq (Qq) broadening cx (Sqj^^Ih ^ and dis- 
sociation, results in the characteristic rapid decrease in 
R^a{pt) hi this part of phase space. In both gold and 
copper reactions at RHIC the suppression -R^^ ~ Raa 
for > 4 GeV and these approach the quenching of 
light hadrons for pr > 10 GeV. It should be noted that 
while for the D mesons we observe a transition from colli- 
sional dissociation to partonic energy loss in the studied 
kinematic domain, for B meson at RHIC the compet- 
ing hadronic processes, Eqs. (|T7)) and ([iSjl . still play the 
dominant role. 

Our theoretical results, presented in Fig. [6l are most 
relevant to the future vertex detector upgrades at RHIC 
that will ensure direct and separate measurements of the 
D and B mesons. These new experimental data will then 
allow to pinpoint the mechanisms of heavy flavor sup- 
pression in the QGP. Indirect studies of D and B meson 
attenuation are currently carried out through the semi- 
leptonic decays of the charm and beauty hadrons. We 
use the PYTHIA event generator [i^ to simulate the 
full kinematics of these Dalitz decays in p+p and A+A 
reactions. The nuclear modification ratio -R^^(pt) of 
inclusive non-photonic electrons is presented in Fig. [7| 
for y/SNN = 200 GeV collisions at RHIC. Predictions 
for central Cu+Cu collisions are also shown for com- 
parison to upcoming STAR data. We remark that the 
Cronin effect included in our study changes the abso- 
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FIG. 8: Suppression of D and B meson production in central 
Pb+Pb collisions at ^/snn = 5500 GeV at the LHC in two 
different scenarios. The top panel shows the quenching of 
heavy hadrons only due to partonic energy loss. The bottom 
panel gives Raa for Ds and Bs with partonic energy loss as 
well as coUisional dissociation of heavy meson. Insert shows 
the corresponding attenuation of non-photonic electrons in a 
limited pr range. 



lute scale of the differential heavy quark cross section at 
intermediate pt- Consequently, to obtain a good descrip- 
tion of the non-photonic -I- e~ quenching results from 
PHENIX M and STAR [H] we use ^ = 3 - 5^ approxi- 
mately 50% larger than in our previous study [3| but still 
compatible with the enhancement of the parton broaden- 
ing that comes from the power-law momentum transfer 
tails of the in-medium Moliere scattering. 

Results for the suppression of open heavy flavor final 
states in central Pb+Pb collisions at ^/sNN = 5500 GcV 
at the LHC in a medium of soft gluon rapidity den- 
sity dNs/dy = 2800 are presented in Fig. E The top 
panel illustrates the effect of partonic energy loss with 
Raa(pt) = Raa{pt) for pT about and above 50 GeV. D 
and B meson quenching approaches the suppression level 
of light hadrons only at very high transverse momentum 
due to the significant gluon fragmentation component to 
pion production in this kinematic region. The bottom 
panel shows the attenuation of open heavy flavor when 
we include collisional dissociation in the medium. The 
most important feature of our results, when compared 
to traditional jet quenching studies, is that at the LHC 



Raa{vt) 



Raa{pt) at a much lower pr — 10 GeV. The 
nuclear modiflcation factor of non-photonic electrons is 
also shown for completeness in the insert of Fig. [S] By 
comparing the top and bottom panels of Fig. [5] we con- 



clude that only experiments at the LHC will have the 
transverse momentum coverage to fully explore heavy 
flavor dynamics in the QGP. We finally remark that the 
small difference in the suppression of D and B mesons 
at Pt ~ 50 GeV is from the interplay of two quenching 
mechanisms - meson dissociation and parton energy loss. 
This difference disappears at pT '^-^ 150 GeV. 



VI. CONCLUSIONS 

Detailed theoretical studies have shown that partonic 
energy loss of heavy quarks in the hot QGP [6|] cannot 
explain the large suppression of non- pho tonic electrons 
observed by PHENIX [H and STAR [l3| collaborations 
in central Au+Au collisions at y'sjvw = 200 GeV per 
nucleon pair. This discrepancy can naturally be resolved 
if collisional dissociation of heavy mesons that tend to 
form inside the QGP is taken into account In this 
paper we reported first results from an approach that at- 
tempts to combine charm and beauty quark quenching 
with D and B meson inelastic breakup processes with 
the goal of describing open heavy flavor production in 
nucleus-nucleus collisions over the full pt range that will 
be experimentally accessible at RHIC and the LHC. The 
treatment of cold nuclear matter effects for massive final 
state mesons was brought on par with that for pions and 
photons psj and the Cronin enhancement at transverse 
momenta ~ 4 GeV was found to be the most important. 
To put studies of heavy meson formation and dissocia- 
tion in the plasma 0, Sj on firmer theoretical ground we 
investigated in the framework of potential models [l3| 
the existence of such bound states in the vicinity of Tc- 
We found that while the temperature at which D and B 
mesons cease to form depends on the details of the in- 
medium quark- antiquark potential and the light quark 
mass, as a rule bound states survive well above Tc. Us- 
ing the light-cone description of hadrons [ll| and the op- 
erator definitions of distribution and decay probabilities 
from factorized perturbativc QCD [s^ we calculated the 
charm and beauty PDFs and FFs in a co-moving plasma. 

In the instant wavefunction approximation, relevant 
to an out-of-equilibrium jet propagation through the 
medium when the timescale for the onset of thermal ef- 
fects exceeds Tform, we evaluated the D and B meson 
cross section suppression from partonic and hadronic in- 
teractions in the QGP. We found that for pT > 4 GcV 
at RHIC and > 10 GeV at the LHC i?f^(pT) ~ 
R^a{pt) and the transverse momentum dependence of 
this attenuation is weak. In our study meson dissocia- 
tion played a dominant role up to pT ^ 10 GeV for open 
charm and up to pt ~ 30 GeV for open beauty hadrons. 
To obtain a good description of the non-photonic elec- 
tron suppression in Au-|-Au collisions at RHIC after the 
inclusion of Cronin enhancement approximately 50% 
larger collisional broadening [2^ that leads to the meson 
breakup was required. Predictions were given for the 
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D and B hadron and the corresponding non-photonic 
electron quenching in Cu+Cu and Pb+Pb colUsions at 
RHIC and the LHC, respectively, that will very soon be 
confronted by data. We conclude by emphasizing that 
while charm quark/meson dynamics in the medium is 
well within the reach of RHIC experiments, for a com- 
prehensive and definitive study which can disentangle the 
medium effects on beauty quarks/mesons not only the en- 
hanced cross sections but also the extended pT reach at 
LHC experiments will play a critical role. 



The solutions to these equations can be written a nice 
form, which we use in the calculation of the heavy quark 
fragmentation functions: 



s) = 



v{p,s) = 



1 



{p-a + m)a 



(A5) 



(A6) 
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APPENDIX A: NOTATION AND 
QUANTIZATION OF QUARK AND GLUON 
FIELDS 

To quantize the theory on the light cone we choose the 



light-cone time as = {l/^/2){x^ 



The light-cone 



"spatial" components are taken as {x , x) where x = 
{l/\/2){x^ ~x^) and x is a two component spatial vector 
with (x)^ = x^ and (x)^ = x^. With this definition, the 
dot product between two four vectors a and b has the 



form a ■ b = a^b 



'6"*" — a • b. We use the notation 



where a stands for the spatial components of a four vector 
a^, and the 3-tuple ~ (a+, a). 

For quarks of every flavor and color there is a corre- 
sponding Dirac field ipcfix). In what follows, we will 
suppress the color index c and the flavor index / for sim- 
plicity. The free field Lagrangian for is given by: 



(Al) 



where ^ = l^d^. 7^ are Dirac matrices satisfying 
{i^tI^} — 27^^"^. We write the Dirac matrices in the 
chiral representation: 



7° = 









^ 








, -a^ j 



The light-cone Dirac matrices are defined by the relations 
7± = (1/V2)(7° ±7^), and 7^ = 7' for i = 1,2. 



The equation of motion for 1/' is simply 

(z^ — m)i]:{x) = . 



(A3) 



We look for plane wave solutions of form u{p) exp{—ip-x) 
and v(p) exp(ip • x) with p'^ > 0. To satisfy the equation 
of motion the four momentum p should be on shell. In 
addition, the Dirac spinors u(p) and u(p) satisfy the fol- 
lowing equations: 

(/i - m)u(p) = , [i) + m)v{p) ^ Q . (A4) 



where is a two dimensional unit vector, which is dif- 
ferent for u and v and depends on the chosen representa- 
tion. For example, based on the choice of and 
u(p, s) and w(p, s) may or may not be helicity eigen- 
states. In our calculations we only use the following 
properties: ^"'""t^"/" = S^s' and Y.s,s' ^sCj ^ss' = -Ux2- 
= (1,1?)'' and = {1,—(t)^ in the instant form, 
where t? are the usual Pauli matrices. 

The general solution to the classical equation of motion 
can then be written as: 

= / -^^iapsuip,s)e-^^-^ + blAP,s)e^'-'') , 

(A7) 

where we will soon identify Ups and bj,^ as creation and 
annihilation operators. 

In light-cone coordinates, not all the components of 
■0 are dynamical, meaning that the time derivative of 
some components do not appear in the Lagrangian. To 
separate these, one defines projection operators A+ and 
A_ as follows. 

A± = ^7^7^ = ^7°7^ - ^7^7" • (A8) 

They are hermitian matrices satisfying A_|_ -I- A_ = 1, 
Aj_ = A±. One defines the "plus" and "minus" compo- 
nents of -0 as the projections "0+ = A+0 and ip- ~ A^-ipj 
respectively. Using the relation 7+A_ = it is easy see 
that '0_ is non-dynamical. Therefore, one can use the 
equations of motion to write it in terms of the dynamical 
field "0+5 which has the following form: 



fori = 1,2, 3. (A2) ip+ix) 



dp^cPp 



(27r)32p 



-{apsU+{p,s)e-'P-^+bl,v+ip,s)e'P-'') 

(A9) 

where u-|_(p, s) = A+u(p, s) and u-|_(p, s) = A^v{p,.s). 
They have norms given by: 

u+{p,s)u+{p,s') = vl{p,s)v+{p,s') = V2p+6ss' ■ 

(AlO) 

The canonical anticommutation relations for the ip-i- 
and ipl fields are: 



{?/>+(.T+,x ,x),tp\_(x+,y ,y)} 
A 



V2 



6{x- -y-)S^^\^-y) 



(All) 
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The anticommutators of two fields or their hermitian con- 
jugates vanish. Writing ip in terms of the a and b opera- 
tors, Eq. (|A9j) . we find that, 

- 6,,,,2p+ {27: fS{p+ - q+)5^'Hp - Q) , 

{bps,bl,} - d,s'2p+{27:)Hip+ - q+)d^'\p - q) , 

(A12) 

and that aU other anticommutator combinations should 
be zero. Therefore a and can be interpreted as anni- 
hilation and creation operators for the fermions, respec- 
tively, and b and 6^ as annihilation and creation operators 
for the antifermions. 

Finally, we describe the quantization of the gauge field. 
The free field Lagrangian is given by: 



C 



(A13) 



where F is the color field strength tensor. The main new 
feature that arises during quantization is that we need 
to fix a gauge. We choose the fight-cone gauge, = 0, 
consistent with the gauge choice for the diagrammatic 
calculations. In this gauge the operator A~ is non dy- 
namical and can be eliminated using the equations of 
motion. 



A- 



(A14) 



The dynamical degrees of freedom are the A fields. 
We define two polarization vectors with perpendicu- 
lar components e(+l) = (l/\/2)( — 1, — i) and e(— 1) = 
(l/-\/2)(l, — «)• The quantized operator A can then be 
written as 



A(x) 



(A15) 



with A = ±1. The canonical commutation relations are 



[5_A'(x+,a;-,x),A^(x+,zj-,y)] 



<5(x--y-)5(2)(x-y) , 



(A16) 



which imply: 

[cpx.cl^,]^6xx'2p+{27rf6ip+-q+)S<^^Hp-q) . (A17) 

Cpx annihilates as giuon and c^^ creates a gluon. The 
four-component gauge field = [^4+ , A~ , A] can be 
written compactly by introducing the four component 
polarization vectors e(p, ±1) = [0, (e(±l) • p)/p+, e(±l)], 
which satisfy p ■ e(p. A) = 0: 

(A18) 



With the creation operators for the particles in hand, 
we can generate the multi-parton Fock states. A single 
quark state can be written as: 

\p^.a)^aUp^)\0) , (A19) 



where the creation operator 



7t 



will carry color. 



spin/helicity/polarization and quark flavor indices a, as 
appropriate. These are normalized as: 

{<t,a'\p^,a) = 5o.a'2p+ {27: fS^'\p^ - t) ■ (A20) 



APPENDIX B: MESON WAVEFUNCTION IN 
INSTANT AND LIGHT-CONE FORMS 

In this section we describe the relation between the 
light-cone meson wavefunction and the meson wavefunc- 
tion in instant form. We do this explicitly for the case of 
the lowest order Fock component that consists of a heavy 
quark Q and a light antiquark q. In the instant form, in 
the hadron rest frame, we write the meson wavefunction 
as: 



91 



d^q2 



1 



-/(gl) 



V3 V2 



(Bl) 



where E refers to the on-shell energy Ei = \/ qf + ml 
and Eh = -y/m^ + is simply mh in the rest frame. 
Here, f{q) represents the momentum space wavefunction 
normalized such that. 



(fig 



1 



(B2) 



In this manuscript we concentrate on distributions 
f{(f') of spherically-symmetric type, relevant to ^5*0 and 
35*1 open heavy flavor states. The particular exam- 
ple for our case will be a harmonic oscillator distribu- 
tion / ^ exp(— (j^ao/2). The invariant mass squared of 
the multi-parton system of no net transverse momentum 
(P = 0) is given by: 



M\q)^2\Y,Pt 



J 




(B3) 



'j-i — "(-j T-Qi ■ We see that is a function of q, 
and the mass of the hadron can be intuitively thought of 
as m\ = M'^{{q)). For the rest of the section we will con- 
centrate on the specific case of a quark-antiquark state 
M2 = {El +E2Y. The form of K'P simpfifies in the 
following two special cases: 



-I- ci(mi, TO2) , if (f '> m. 



2 



(B4) 



(tTti+m2)^ 
(Tnim2) 



+ C2(TOi,m2), if < . (B5) 
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In all cases the same functional dependence on k = q 
and X is obtained via ~ and the numerical coef- 
ficients Ci{mj) are absorbed in the proper normalization 
of hadron state. This yields: 



if) ^ cxp 



+ (1 — x)ml + xm\ 
2A2x(l - x) 



(B6) 



In two limiting cases the coefficient that relates A in 
Eq. (jB6P to oo (the width of f{(f)) can be analytically 
estimated, as discussed above. Since the form of -0 is the 
same in the two extreme cases, we find it more conve- 
nient to treat A as a constant which can be determined 
through the requirement that we obtain the same expec- 
tation value for the transverse momentum from the in- 
stant and light-cone wavefunction forms, i.e. (k^) = (q^). 
For the example of typical open heavy flavor wavcfunc- 
tions Q the numerical values arc within 30% of the values 
for the two limiting cases (|B41 IB5|) . 

To establish the correspondence between the instant- 
form wavefunction and the light-cone wavefunction, we 
now convert the spatial momentum integration variables 
q to ~ (q"*" , q) . The Jacobian for the transformation 
is q'^ /E{q). To recover the (anti)commutation relations 
outlined in Appendix lAl we find that no additional factors 
arise for the creation and annihilation operators. The 
proper normalization of states oc q+ is ensured when we 
observe that (5(q^ — ) = ((7+ /£■((/)) (5 ((7+ — g+ ). Similar 
relations can be established for antiquarks and gluons. 
One can also rewrite 5[ql + f/f — P'^) {■\/2) 5{qf + 
q2 — P'^) in the boosted frame where we will match the 
instant and light-cone forms of the state. P+ is the plus 
component of the light-cone momentum of the meson. 
Note that in the meson rest frame = 0. Making these 
variable changes and recalling that q = k leaves us with 
the same meson state: 



\P' 



^y(2^^/25; 



si^x,-! 



xf{fix,,k,))al^^'^iq-{ + )bl^'^^Hq.2 + m . (B7) 
We complete the correspondence by observing: 




hm 2VE„J^f{q'{xM) 



P+V2f{x„k,) = ij{xM . (B8) 



The prefactor of V P+ in front of the rest frame wave- 
function f{k) in the relation between ip{xi,]ii) and f{k) 
(Eq. IBSP can be physically understood as follows. A 
boost from the rest frame to a frame with large momen- 
tum P+ leaves the transverse shape of the wavefunction 
unchanged. However, in the large light-cone momentum 
direction the wavefunction is stretched in momentum 



space by a factor V P+. Note that Eq. 
with Eq. (IB2| yield J d'^kdx\^l!{x,k)\'^ = 
ployed in this paper. 



(IBSP combined 
2(27r)^ as em- 



APPENDIX C: MEDIUM CONTRIBUTION TO 
HEAVY QUARK FRAGMENTATION 
FUNCTIONS 

In Sections IIIII and [TVl we calculated the fragmentation 
functions for c and b quarks into D and P-mesons re- 
spectively, from the definition given in [s^ , Eq. . In 
the adiabatic picture, where the meson wavefunction has 
sufficient time to equilibrate with the medium, we take 
the modification of D and B meson wavefunctions in the 
medium into account by using the thermal wavefunction 
calculated from potential models in Section IIVI 

In this equilibrium situation, one can attempt to go 
a step further by estimating additional corrections that 
arise from thermal parton pick-up as follows. Instead of 
taking the expectation value in vacuum, we can define 
this thermal fragmentation as: 



D 



1 7+ 

^^'^color^^'^D^rac^ (CI) 



Y.ra eM-l3Em){m\^{y-,0)al{P+)ah{P+mO, 0)|m) 
E™exp(-/3i?„) 

We take Eq. (jCip to be the minimal generalization of the 
definition of fragmentation in perturbative QCD, with 
z < 1. In Eq. (|C1|) . Em is the energy of the state |to), 
13 is the inverse temperature, and exp(— /3Pm) gives the 
thermal weight of the state m. For convenience of no- 
tation, we are writing the spectrum of states as being 
discrete, with the states normalized as (m|n) = 5mn- 
One can immediately see that in the case that the op- 
erator does not contract with the state to, we get back 
the definition in vacuum because the sum of exponentials 
cancels out between the numerator and the denominator. 
For weakly interacting quasi-particlcs, is the sum of 
the energies of the quasi-particlcs in the state m. In this 
approximation, the denominator 



Z = Y, cxp(-/3S„,) = n,, (exp(/3e,) ± l) 



(C2) 



where the -I- is for Fermions and — for the bosons. When 
the operator contracts with one quasiparticle with energy 
e from the medium (equivalent to saying that the meson 
picks up a thermal parton of energy e), all the other 
terms in the product ('Eq. lC2[) cancel out as before, except 
the one corresponding to energy e. Hence, one gets an 
additional factor of g{e) = l/(exp(/3£) ± l) depending on 
whether the parton is fermion (+) or a boson (— ). 

At the same order in as vacuum fragmentation the 
following Feynman diagram is the most important addi- 
tional correction to Dh/Q{z) function. Fig. [HI We call this 
correction P/j/q^^-* (z). This diagram can be interpreted 
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FIG. 9: Thermal correction to the fragmentation function 
that arises from an in- medium quark pick-up. 



as follows. The heavy quark picks up a light thermal 
antiquark from the medium and forms a meson. The 
evaluation of this diagram is very similar to Eq. (j26p and 
yields: 



dxid'^kiip{xi,'ki) d.T2d2k2?/'*(a;2,k2) 



V2 V2 



(27r)32Va:;i(l-a;i) (2^)32^X2(1 - .T2) 



ds 9 



1 - z 



Tr 



7 • Pa - "1( 
A 



-7^Usi {pi)Vs' {p'lh^i-J ■ Pq + TOg)7'^ 
) 

■np^(pQ - Pg)ncrA(P6 - Pq) 



^ -Pb- m, 



9iP°q) 



1 



Tr[7+(7 -p)] ' 



(C3) 



where Pa = Pb = Pq and pq is the momentum of 
the light antiquark that is picked up from the medium. 
ff(Po) = 3„u is the thermal probability to find a quark 



^1+1 



with energy p^^ in the medium. For future use we define 
the hadron momentum ph = pi + p'l = P2 + p'l- The ex- 
pression, Eq. (|C3p . differs from D}^/q{z) only in the extra 
factor (?(Pg). To get a rough estimate of the correction, 
we replace p° by its typical value p° and obtain. 



i?W(z)~i?(z).g(p°) 



(C4) 



^From the Feynman diagram Fig. ^ we see that Pq = 
Pa Ph- For a typical sample of the kinematics, we con- 
sider the following. By choosing an appropriate direction 
of the axes, the transverse momentum of the hadron can 
be taken to be zero. Then, 



Pa 



Ph 



Pq 



P 



Pa, Pq 



2zp+' 
{l/z-l)P 



(C5) 



2F+(l/z- 1 



-,P<j 



where Pg is the transverse momentum of the light anti- 



quark and is related to the off-shcUness [pj^ - 



of the 



initial heavy quark. We obtain straightforwardly: 



fq = 



1 r 
71 



P^{--i) 

z 



2P+(l/z- 1) 



(C6) 



>^P-(1/-1) 



where is the light quark mass in the thermal medium. 
To obtain a conservative estimate for the correction, we 
take mq to be zero. Using typical value of z for D and 
_B-mesons, we find that for a 4 GeV D-meson moving 
through the medium at a temperature 250 MeV, this 
gives g{fPq) less than 3%, while for a 4 GeV i?-meson 
moving through the same medium we find corrections 
less than 15%. 



APPENDIX D: ANALYTIC MODEL OF HEAVY 
QUARK AND HEAVY MESON MULTIPLICITIES 

It is useful to find an approximate analytic solution to 
the system of ordinary differential equations, Eqs. ((TT]) 
and ([35]), for the purpose of understanding the inter- 
play between heavy quark fragmentation and meson dis- 
sociation and evaluating the fraction of the time dur- 
ing which the charm is in a partonic state. This will 
help us determine the quenched initial conditions for 
our numerical results. We will take the pT-dependent 
quark and hadron multiplicities, f'^{pT,t) = Aq/v^ 
and f^{pT,t) = Ah/px" to be of power-law type [2^. 
We can then express: 



dx — 

X' 



4Q/H{x)f"{pT/x) = f"{pT/cH) , 



and 



dz -DH/Q{z)f'^{pT/z) ^ /"^(pt/cq) 



For any finite pr range the coefficients uh, "q, ch, cq 
can be determined numerically. In order to obtain a sim- 
ple analytic solution we must also average the formation 
and dissociation rates in the interval (0, Lqgp)- This 
yields: 

OtJ^[PT,t) - 7 TT H -, -, 7^, (Ul) 

(Tform (PTj) (Tdiss(PT)) 

dtf {pT,t) = — y- + ^ — (D2) 

A general solution ansatz oc ^ Aie""'* leads to: 



1 



1 



(Tform) (Tdiss) 



(D3) 




'^form'^diss 



-(1 
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rvalues From the initial conditions, Eq. (|49p. and 
Eqs. (|D1[) and (|D2p we can also determine the time 
derivatives of the heavy quark and hadron distributions 
at i = 0. Solving for the cocfhcicnts in the linear super- 
position of independent solutions we find: 



f"iPT,t) 



V form 



{riormiPT)){ri - r2) 



(D4) 



/'5(PT,0) 



and, 



1 



''"l ^2 ^ ('^form 



r-2 



ri - r2 MTtorm(PT)) 



(D5) 



We can use Eqs. (|D4p and (jPSP not only to study 
the mix of heavy quarks and hadrons as a function of 



('''form(PT)) , (Tdisslpr)) but to also evaluate the fraction 
of the time the charm or bottom arc in a partonic state: 



t 



partonic 



(pt) 



itotal 



(D6) 
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